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Abstract 

A new k — e turbulence model that accounts for viscous and wall effects is 
presented. The proposed formulation does not contain the local wall distance 
thereby making very simple the application to complex geometries. The for- 
mulation is based on an existing k - e model that proved to fit very well with 
the results of direct numerical simulation. The new form is compared with nine 
different two-equation models and with direct numerical simulation for a fully 
developed channel flow at Re = 3300. The simple flow configuration allows a 
comparison free from numerical inaccuracies. The computed results prove that 
few of the considered forms exhibit a satisfactory agreement with the channel 
flow data. The new model shows an improvement with respect to the existing 
formulations. 


1 Introduction 

The Reynolds averaging of the Navier-Stokes equations ( N — S) allows solving very 
complex flow configurations with a manageable computer effort [1, 2]. A possible 
alternative approach is the so called Direct Numerical Simulation[3 ] , (DNS), that, 
because of the current memory and speed limitations of supercomputers, is restricted 
to simple flows and low Reynolds numbers (Re). So, the solution of the Reynolds 
averaged N — S equations still represents the main possibility for flow simulation 
in engineering applications. The well known drawback of the Reynolds averaging 
technique lies in the introduction of further unknowns, the so called Reynolds Stresses , 
stemming from averaging of the non-linear convective terms. 

The Reynolds stress components are normally related to the mean flow quantities 
through a set of additional equations that represent the turbulence model. It is 
possible to carry out a further simplification following the Boussinesq assumption in 
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which the Reynolds stress components are expressed in terms of the mean velocity 
gradients as follows: 


(dUj dUA 2 

= ■" 1 fc- + Wj ) - 3 M " 


in which k = is the turbulent kinetic energy, u; is the fluctuating velocity, 

V{ is the mean velocity, vt is the turbulent viscosity that represents the proportionality 
parameter between the Reynolds stress and the mean strain rate, and the overbar 
stands for the average operator. The turbulent viscosity v% is generally related to the 
local velocity and length scale of turbulence. 

Accordingly, it is possible to define two main model categories, with decreasing 
complexity, that have an acceptable degree of generality: 


• Full Reynolds Stress models. This approach does not use the Boussinesq 
assumption and every component of the Reynolds stress tensor is computed 
using a transport equation. Due to the high non-linearity usually connected 
to the source terms of the differential system, this class of models requires 
numerical methods which can effectively treat the stiffness introduced by the 
source terms. 

• Two-Equation models. In these lower order models the Reynolds stress com- 
ponents are computed taking advantage of the Boussinesq assumption so that 
the degree of complexity of the formulation is considerably reduced. The turbu- 
lent viscosity is normally computed using two transported variables related to 
the velocity and length scale of turbulence. The most popular choice of variables 
is the turbulent kinetic energy k and the turbulence dissipation rate e defined 
as e = vuijuij. 

In their early forms, both classes of models treat the solid boundary and molecular 
viscosity effects by using the wall function approach[4]. In this method the viscous and 
buffer layers in the wall proximity are modelled by an algebraic expression that relates 
the turbulent quantities to the shear velocity under the hypothesis of local equilibrium 
turbulence. While the assumption is valid for fully developed flow conditions, it is 
the well known cause of inaccuracies in case of strong streamline curvature, adverse 
pressure gradients or in stagnation regions that are often encountered in practical flow 
configurations. The full Reynolds stress models have been applied to the study of an 
impinging jet in a squared duct [5] with satisfactory agreement with experiments, 
whereas two equation models have been used for the simulation of a much wider class 
of flows[l]. Unfortunately, most of the applications have been carried out by the 
high Re formulations that still retain the wall function approach in presence of solid 
boundaries to avoid the difficulties connected to the simulation of the viscous and 
buffer layers. 

More recently, with the introduction of new forms, called low Re y (XiZ), of both the 
full Reynolds stress and two equation models, it is possible to account for the molecu- 
lar viscosity and the wall effect [6, 7]. The LR forms differ from the standard high Re 
forms insofar as they include further terms to model the influence of molecular vis- 
cosity that is not negligible in wall regions, the normal velocity fluctuation damping 
exherted by a solid boundary, and the presence of a nonisotropic contribution to the 
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dissipation rate of turbulence that becomes dominant in the viscous layer. The LR 
models have been designed to maintain the high Re formulation in the log-law region 
and further tuned to fit measurements in the viscous and buffer layers. Because of the 
difficulties in measuring the turbulent stresses in the wall proximity the models have 
been made to reproduce measurements with a very wide uncertainty range, thereby 
making turbulence model refinement very unreliable. The recent availability of the 
results of direct numerical simulation of turbulent flows, even if restricted to very 
simple flow configurations and quite low Reynolds numbers, helps in understanding 
the limiting behavior at boundaries of the turbulence quantities and provides a valu- 
able data base for model testing and refinement. With this in mind, it is possible to 
improve simple models making them consistent with the abovementioned data set. 

In this paper the attention will be focused on lower order two-equation models in 
their LR form. From a computational point of view, the modifications introduced to 
model the wall effect often cause serious numerical problems so that the results of the 
computations may be heavily affected by the adopted solver. A comparison of the 
predictions given by the existing k — eLR models obtained with a single solver able to 
cope with all the formulations is definitely desirable and wets one of the targets of the 
present investigation, together with the development of a new general k — e model. 
This also allowed investigating the numerical behavior of the models. 


2 k — e LR models 


In the current literature it is possible to find several wall effects modifications to the 
standard high Re k — c model. Among the various LR formulations nine of them were 
considered for the comparison with the present model (pr) proposed in this paper: 
Chien[8], (c/i), Jones and Launder[9], (j/), Nagano and Hishida[10], (n/i), Coakley[ll], 
(co), Speziale et al.[12], (sp), Kim[13], ( ki) } Rodi[l4], (ro), Lam and Bremhorst 15 , 
(/&), Shih[16], (sh). In parentheses are reported the letter codes with which the models 
will be referred to. The selection was performed on the basis of the results of Patel 
et al. [6] with the addition of more recent models. 

The LR models may be expressed in a standard form in terms of the following 
nondimensional variables (the overbar stands for a nondimensional quantity): 


k = 


k 

U2 


€ L V t 

? = U3 *'*“7 



Xi = 


ft 

L 


in which U is a typical velocity, L is a typical length, and v is the molecular viscosity. 
Accordingly, the flow Reynolds number is defined as: 


Re — 


L U 

v 


Hereinafter the variables will appear only in their nondimensional form so that the 
overbar will be dropped for simplicity. The model transport equations may be written 
as follows: 


dk , 
Si + 


^r = *(3§;( 1 + ft)S) + p -' + c + n 

If + = -k (A 0 + £) &) + «./. iP - 


d + E 


( 2 ) 


3 


in which 


dUj dU % 
dx{ dxj 


P - 




is the production rate, and 


k 2 

Vt = Re — 


( 3 ) 


is the Kolmogorov expression for the turbulent viscosity. The extra terms D, E, II 
represent the corrections to the standard high Re formulation that, together with the 
damping functions f\ y / 2 , / M , allow balancing the transport equations (2) in the wall 
region. 

The term D represents the extra destruction rate in the wall region. The dissipa- 
tion rate appearing in (2) may represent the total dissipation rate, e, or the so called 
isotropic dissipation rate, e, that are related via the following expression: 


€ = € + D (4) 

in which D is the dissipation budget in the wall region. Evidently, e = l far from solid 
boundaries where D — ► 0. A similar correction is generally made in the e transport 
equation. The extra term II is the pressure transport term for the trace of the 
Reynolds stress tensor. Shih[16], analysing the DNS computations, observed that II 
is o(y) in the wall region, but is generally neglected in all the existing two-equation 
models although the turbulent transport is o(y 3 ) and remains in the model equation. 
This term is modeled as a turbulent kinetic energy diffusion. 

In all the proposed formulations, the damping functions /i, /^, are related to the 
Van Driest damping expression for the mixing length, L : 

L = y(l-e~ Ay+ ) 

in which A = 26 is an experimental constant, and y + is a nondimensional wall distance 
defined as 

y + — u T y Re (5) 

The selected transported variables are k and e with two exceptions: 

• co. The transported variables are q = y/k and u> = p According to Coakleyfll], 
this choice ensures better numerical properties, as explained in ref. [17]. 

• sp. While maintaining k as the first variable, e is replaced by r = | in order to 
remove the inaccuracies connected with the solid wall boundary condition for 
the dissipation rate, since r vanishes automatically to zero. 

The to two-equation and ki four-equation models are two-layer models in which 
an inner layer is defined close to the wall ( y + < 100) where, while using a transport 
equation to compute the turbulent kinetic energy, the dissipation rate is computed 
by algebraic expressions in function of the wall distance and k . This choice allows 
removing most of the numerical problems since the cumbersome solution of the dissi- 
pation rate equation together with the implementation of a proper boundary condition 
are skept in the wall region. Moreover, no extra terms are needed in the transport 
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model 

Jp 

/ 1 

/ 2 

ch 

1 — exp{— .0115?/ + ) 

1.0 

1 - 0.22e*p(-§) 



1.0 

1 - 0.3 exp(—R 2 ) 

nh 

0 - ex P(-^s)) 2 

1.0 

1 - 0.3ea;p(-.ft 2 ) 

CO 

1 - exp(- .0065i? y ) 

1.0 

1.0 

sp 

o + - ianh (^)) 

1.0 

1 - 0.22ezp(-§) 

ki 

1 - exp(—.025\/Ri — 10~ 5 R?) 

1.0 

1.0 

ro 

1 — exp{— .0198^) 

1.0 

1.0 

lb 

(l-ea:p(-.0165^)) 2 -(l + ^r) 

(1+f) 2 

1 - exp(-Rf) 

sh 

1 - exp{-Y?i=\ a iy + ') 

1.0 

1-0.22 exp(-§) 

pr 

. exp(— .0004erp(1.2/?^ 4 )) 

1 erp(-.0004) 

1.0 

1 - 0.22exp(--§) 


Table 1: Damping functions 

in the LR forms 
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E 
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0 

~R~eJ 

7 ^exp{-. 5 y + ) 

ji 

0 

_j_ (d^k ) 2 
Re \ dy ) 

24 / a 2 t/\ 2 
R? \ W ) 

nh 

0 

J2_ ( d\Lk\ ^ 

Re \ dy ) 

-^(J - U) (0) 
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0 

— 
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0 

- 
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0 

— 

— 

ib 

0 



sh 


0.05 /it jL 1 

0 

ill— (d?U ) 2 

Re 2 V dy 2 ) 


/^l-erp(-t/+)) a k K ^\ j 

pr 


It 

0 

^ (dlu ) 2 

Re 2 \ dy 2 / 


Table 2: Extra terms in the LR forms 
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equations since the wall effects are supposed to be effective only within the range of 
application of the algebraic expression for e. 

The damping function expressions selected in the various formulations may be 
found in table 1, whereas table 2 gives all the selected forms for the extra terms D, 
E , II. Most of the models have a set of damping functions based on the physical, y, 
or nondimension al, distance, and on the turbulent Reynolds number based on y, 

R y = y Re , or R t = The exponential function introduced in (3) damps 

the turbulent viscosity to zero in the viscous and buffer layers. It is interesting to 
observe that nearly all the models retain the exponential function for the decay of 
dissipation rate f 2 proposed by Hanjalic and Launder[18], according to which c 2 has 
a finite value at the wall. The sp formulation needs an additional function to further 
damp the value of c 2 in the wall region (see ref.[12]). The empirical constant ci is 
left unchanged all the way down to the wall in all the formulations with the only 
exception of lb where, since no extra terms are added in both k and c, the production 
of dissipation rate must be damped to balance c at the wall. 

For any further detail about the models, the reader can refer to the original 
bibliography. 


3 A new form independent of wall distance 

The results of DNS have shown the limiting forms of the turbulent quantities in the 
wall region. A correct formulation should ensure that 

U = o{y ) k = o(y 2 ) e = o(l) 

l = o(y 2 ) fi t = o(y 3 ) uv = o(y 3 ) 

together with a negative slope of the dissipation rate in the wall region. The two- 

equation model proposed by Shih[16] matches the flow variables limiting forms re- 
markably well, but in its original formulation suffers from numerical problems stem- 
ming from the selected form of both the dissipation rate decay term and turbulent 
viscosity formula. In fact, the dissipation decay rate is 


r & 

- C2 H 

that is o(l) at the wall, and the turbulent viscosity is computed as 

k 2 

v t ~ Cn fa Re — 

€ 


( 6 ) 


( 7 ) 


in which the damping function is based on a fourth order polynomial of t/ + designed 
to match the DNS data. The isotropic dissipation rate is computed as 


e = e — 


2 k Re 


( 8 ) 
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This expression ensures that c vanishes at the wall like y 2 , provided that a proper 
boundary condition makes € equal to the negative term in (8). Even in the com- 
putation of simple fully developed channel flows some convergence problems were 
encountered, which were caused by the decay term given by equation (6). When 
computing the isotropic dissipation rate with equation (8) the ratio of the first order 
derivative of k and k in the wall region often gives overshoots that cause negative 
dissipation rates and, through equation (7), negative v t . 

Alternatively l may be defined as 


€ = cf € (9) 

where we have introduced a function f c which is a function of the turbulent Reynolds 
number, Rt = 

ft = f(Rt) = 1 - exp(-s/Ri) 

At the wall, R t = o(y 4 ) so that the function f e is o(j/ 2 ). Equation (9) prevents any 

negative dissipation rate, makes e = o(y 2 ) at the wall, and is balanced by molecular 
diffusion and the extra production rate E. The f t function is designed to give l « e 
for ?/ + > 6. 

Table 1 shows that all the LR forms except the jl model use damping functions 
for v t based on y or y + , while / 2 is based on R t . The use of as the exponent of the 
damping functions gives unphysical results in the case of stagnation points. In fact, 
the definition of the friction velocity u T appearing in equation (5) is: 


u T = 




uall 

P 


T~wall — Plam 1 a ) 

V on J waii 

where n is the direction normal to the wall. At separation and reattachment points 



0 


and — ► 0 regardless of the value of y. This implies that is zero all along the 

flow section making the viscous layer thickness unphysically unbounded. There are 
several known tricks to overcome this problem, like relating u T to the k peak in the 
viscous layer via the wall function for the turbulent kinetic energy: 


u T 



maxj 


or replacing the velocity gradient by the 
cross flow direction: 



maximum value of the vorticity u in 

~| w ma x | 


the 


While these, and other, tricks remove the singularity, they are somewhat arbitrary 
and represent sources of inaccuracy. Moreover,, a general turbulence model should 
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not need information about the flow domain geometry, like the wall distance, the 
determination of which is not straightforward in complex geometries. Unfortunately, 
replacement of y + is not an easy task. In fact, while becomes constant for y + > 
5 — 10, /p, related to the wall effects, vanishes gradually in the range 0 < y + < 100. 
Figure (1) shows the R t profile as a function of the turbulent Reynolds number 
increases very steeply for 0 < y + < 50, reaches a peak at y + « 80 and decreases for 
y + > 100. When keeping as a benchmark the f ^ = f(y + ) given in [16], which was 
carefully designed to give the best fit with DNS , it has been impossible to replace y + 
by Rt even using complicated hyperbolic expressions because the turbulent Reynolds 
number does not behave monotonically. After intense numerical testing we found 
convenient to introduce a new damping function independent of y and based on 
length scales. Figure (2) shows that the turbulence length scale, defined as 


L = 



€ 



€ 


has a monotonic behavior so that it may be conveniently used to replace y + . A new 
nondimension al parameter, Ri, based on L is introduced: 



( 10 ) 


in which 




\U\ 


| U | is the amplitude of the relative mean velocity in a frame fixed to the solid 
boundary. This choice ensures the Galileian invariance of the model equations. Ri 
represents the ratio of the turbulence length scale L to the viscous length scale L v . 
This ratio approaches zero at the wall because L — > 0, L u — » oo, and reaches asymp- 
totically a maximum in the log-law region. In the wall region the following limiting 
forms hold: 


L = o(y 3 ) 
U = o(y) 


Rl = o(y 4 ) 


1 

Accordingly, the exponent of the damping expression must be R A L to give = o(y) 
at the wall. The new damping function based on (10) has the following expression: 


/m = 1 - exp (- c ) • eX P ( _c ^ eX P ( c ^ 2 R L 4 )) ( n ) 

in which c^i and 2 are two empirical constants. The f ^ function given by Shih[l6] 
versus R l is shown in figure (3). The shape of the exponent of f li (= —ln{ 1 - / M )) 
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closely resembles an exponential: this was the reason why a double exponential func- 
tion was selected. 

Figure (4) shows that the differences between the original damping function[16] 
and the new form are limited to the viscous layer (j/ + < 10), while in the buffer and 
log-law region the new form independent of the wall distance quite closely resembles 
the formulation based on the nondimensional wall distance. 

With the substitution of equation(8) by equation (9) and the replacement of the 
damping function for the turbulent viscosity by equation(ll), a new model tuning 
was necessary. In the present formulation II has the form given in table 2, c\ = 1.45, 
c 2 = 2.0, <Jk = &c = 1-3, c M i = 0.4 * 10” 3 , Cp 2 = 1.2. The boundary condition for the 
dissipation rate is 


^ wall — 



2 k Re 


( 12 ) 


4 Fully developed turbulent channel flow 

Patel et al. [6], while comparing some LR forms with experiments, stressed the data 
unreliability in the wall proximity in the range 0 < < 100. In this region the 

scatter in the turbulent kinetic energy k is of the order of ±30%, and comparable 
inaccuracies may be expected for the turbulent shear stress uv, defined in equation 
(1), and the dissipation rate 6. Because of these uncertainties it was decided to use the 
direct numerical simulation results of a fully developed channel flow by Mansour et 
al.[3] at Re — 3300 based on the maximum velocity, or J?e r = 180 based on the friction 
velocity. All the derivatives in the flow direction are zero, with the exception of the 
streamwise pressure gradient that is constant. This allowed studying the problem 
by a one dimensional grid with 81 points in the cross flow direction and a geometric 
stretching ratio of 1.05. With this grid it was possible to place the first grid point away 
from the wall at j/ + « 0.15. The transport equations are nondimensionalized with 
respect to half channel height and the friction velocity u r , while the turbulent viscosity 
is made nondimensional with respect to the laminar viscosity v . The domain and the 
solution are symmetric with respect to the centerline where the following boundary 
conditions were imposed: 

dV_ _ dk _ <h 

dy dy dy 

Under these assumptions continuity is automatically satisfied and the momentum 
equation in the main flow direction, x, is simplified to: 



in which y represents the cross flow direction and U is the flow velocity. The transport 
equations for k and e become: 

B^A ( 1 + £)f)+ i> - E + fl + n = ° 
t fc (£(» + £) -<*/*£ + £ = o 
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in which the one dimensional production rate is expressed by: 



The differential system is solved by a decoupled approach based on a simplified 
approximate factorization algorithm described in ref. [2, 17]. The comparison of the 
various models is done by plotting the mean velocity £7, the turbulent kinetic energy k y 
the dissipation rate c, the turbulent shear stress u\ 7, the turbulent viscosity and the 

production rate P. The plots gather models that solve for the isotropic dissipation 
rate e, (ch y jl, nh) y models that use a different set of variables (co, sp), two-layer 
models (&i,ro), and models that solve for the total dissipation rate ( lb y sh,pr ). The 
sh results reported here have been obtained with a slightly modified version of the 
model in which the original expression for l given in equation (8) was replaced by 
equation (9). 

Figure 5 shows the computed velocity profiles compared with DNS , whereas fig- 
ure 6 shows the comparison of the turbulent viscosity profiles. Among the models 
that solve for the isotropic dissipation rate nh is the one in better agreement with 
DNS . This is strictly related to the good reproduction of the turbulent viscosities 
in the viscous and buffer layers. This produces the correct velocity gradient in the 
buffer layer together with a slight overprediction of the centerline velocity. The good 
agreement is obtained by nh despite the fact that, at the wall, u t = o(y 4 ) because 

= o(y 2 ). jl gives a strong underestimation of the centerline velocity. This is 
caused by the underestimation of the turbulent viscosity in the buffer layer that does 
not ensure enough cross flow momentum diffusion. The reason for this is the damping 
function for the turbulent viscosity / ;i . This function, designed to work for high Re 
flows, in this case does not reach unity for y + > 100 and is still approximately 0.65 
on the centerline. This exceedingly smooth behavior indicates that the selected form 
of the damping function for the turbulent viscosity, while being based on R t and not 
on is not general, sp shows an excellent agreement with experiments (figure 5). 
This is strictly linked with the good reproduction of the turbulent viscosity profile 
(figure 6). In fact, and experiments[19] agree in locating the peak of turbu- 

lent viscosity at 100 < y + < 150 and not on the centerline: this feature is correctly 
reproduced by sp , while is completely missed by co in which the overestimation of 
the turbulent viscosity causes a very smooth and overpredicted velocity. For the two- 
layer models the influence of the matching point was investigated. For ro the two 
computations carried out matching at y+ = 50 and j/ + = 100 showed a considerable 
change in the turbulent viscosity, but negligible changes in all the other quantities 
including the velocity profile. The four-equation, two-layer ki model retains only 
in the inner layer, while it assumes that in the outer layer = 1. The selected 
form of the damping function seems to suffer the same problems encountered with jl 
insofar as at matching (y+ = 100) the damping function is still 0.7. This is shown in 
figure 6 where ki refers to the original formulation by Kim[13], while in ki(fmu) the 
damping function was retained in the outer layer. As it was found for ro, the velocity 
profile is not affected by changes in turbulent viscosity taking place at > 100. The 
three models that solve for e give the best overall agreement with the DNS velocity 
profile. The present formulation, pr, fits with DNS in the buffer layer, and predicts 
the correct velocity at the centerline. The improvement with respect to j7, which is 
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the only other model that does not use the wall distance, is quite remarkable. The 
v t profile exhibits the correct trend, although the turbulent viscosity value is slightly 
overpredicted on the centerline. 

The turbulent kinetic energy profiles, reported in figure 7, show a quite wide scat- 
ter in the computational results. For all the models k oc y 2 at the wall. Nevertheless 
the k peak is either misplaced or underpredicted. An accurate choice of the constant 
in the pressure transport term together with a correct dissipation rate budget in the 
wall region allows the sh and pr formulations to fit both the location and the mag- 
nitude of the turbulent kinetic energy maximum value, although the k level in the 
log-law region, which influences the turbulent viscosity distribution, is overestimated. 
Both the two-layer forms do not seem to match the data not only in the region of 
application of the algebraic expression for e, but also in the outer layer, sp , while 
giving an exceedingly low k peak, predicts the correct level of turbulent kinetic energy 
in the log-law region. 

The turbulent shear stress profiles are given in figure 8. This quantity is critical 
for the correct reproduction of the effect of turbulence on the mean flow field. Not 
all the models satisfy the relation Uv oc y 3 at the wall. In particular, nh y co, and lb 
predict uv oc y 4 . Since all the models ensure a the disagreement is caused 
by an improper slope of the turbulent viscosity induced by the damping functions. 
When both k and e are o(y 2 ), f ^ must be o(y): the latter condition is not satisfied 
by nh and co, in which = o(y 2 ), while for lb the inaccuracy comes from = o(l) 
together with the use of the total dissipation rate e = o(l) to compute the turbulent 
viscosity, sp and pr show the be st fit with experiments. Nevertheless the present 
formulation manages to give the proper uv shape. This has been accomplished by 
optimization of the c^\ and c ^ constants. 

Both the two-layer models show a kink in the turbulent shear stress profile located 
at the interface between the inner and the outer layers. The kink is caused by an 
imperfect match between the inner layer, where an algebraic expression for e holds, 
and the outer layer where the standard high Re transport equation for € is imple- 
mented. The two layers are normally interfaced in such a way that the value of the 
dissipation rate obtained by the algebraic expression equals the one obtained by the 
transport equation. At the interface point the two layers may have the same e value, 
but evidently have slightly different slopes that induce the kink. Different damping 
functions, choice of constants, and matching criteria did not solve the problem, so 
that it is reasonable to conclude that, up to now, it has been not possible to ensure 
a Cl continuity for e at the matching point, and it is actually questionable whether 
it is possible. In the elliptic flow solver based on tl[ 2] the layers were matched where 
vtjv « 20 — 36. In the present low Re flow, this condition is never satisfied so that 
the layers are matched according to the nondimensional wall distance y 4 *. Part of 
the inaccuracies shown by tl may be attributed to this technique which brings the 
matching point too close to the wall. 

DNS shows that the e peak is located at the wall and the slope of the dissipation 
rate is negative in the viscous layer. These two results contradict what have been 
found so far by experiments. The reason for this disagreement is largely attributed to 
difficulties in measuring turbulent velocity gradients close to the wall with the result 
of an extremely large experimental uncertainty. All the older formulations have been 
tuned to have a positive or zero e slope at the wall. The only models that took 
advantage of DNS are sp , sh , and the present formulation: this should be kept in 
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mind when looking at the results. Figure 9 shows the dissipation rate profiles for 
the ten models. For the formulations that solve for e the plots refer to the isotropic 
dissipation rate without the extra dissipation D (equation 4). This term, while being 
dominant in the viscous layer, gives values of e at the wall that are at most half of 
what DNS prescribes, together with a positive slope. The local peak of dissipation 
rate away from the wall takes place at y + zt 12 and is not influenced by the extra 
term D which decays very quickly and is not negligible in the range 0 < y+ < 10. In 
all the models, including sp , sh and pr there is a marked tendency to overestimate 
this peak, with the exception of ch which fits very well with DNS in the range 
10 < y + < 180, probably because of the selected form of the extra production term 
E. The only two models that exhibit the correct dissipation rate slope at the wall are 
sh and the present formulation pr based on the sh model. The correction adopted 
to the isotropic dissipation rate formula in sh , based on an exponential damping 
function, while removing numerical troubles, produces an overestimation of e at the 
wall of the order of 2.5 times the value given by DNS . Although e is o(y 2 ) at the 
wall, it seems that it has been not possible to obtain the proper balance between the 
molecular diffusion of dissipation and decay rate. The problem may be generated 
by the difficulties in the implementation of a proper boundary condition for e. In 
fact, the dissipation rate must reach a finite value at the wall which comes from the 
balance of the abovementioned terms. The errors in computing first or second order 
derivatives at the wall may be the primary__s_ource of inaccuracies in the treatment of 
solid boundaries. The dissipation rate boundary condition given in equation (12) was 
retained for both sh and pr. The overestimation of the wall value of e is reduced to 
40% by the present formulation. It has to be observed that sh and pr are the only 
models, among those considered in this comparison, that show a negative e slope at 
the wall. 

Figure 10 proves the close link between the velocity profiles, together with the 
turbulent shear stress, and the production rate P . The plots show that the models 
locating P ma x at y + a 10 - 11 reproduce very well the velocity profile. If P max , re- 
gardless of its value, is located at y+ < 10 the velocity module on the centerline will be 
underestimated, while shifting P ma x further out generally produces an overestimation. 

5 Conclusions 

A new k — € LR model is presented. The formulation was obtained starting from 
the LR form proposed by Shih[16]. The comparison of the new formulation together 
with nine other LR forms shows an overall improvement in the fit with the direct 
numerical simulation data. The lower order model generality is greatly improved by 
introducing a new damping function for the turbulent viscosity which does not need 
information about the flow geometry. This is done by a new form of based on the 
ratio of the turbulent and viscous lenght scales Rl that removes the uncertainties 
and inaccuracies stemming from the need to define a geometrical distance from solid 
boundaries. It is necessary to test the new formulation independent of y+ for a wide 
range of Reynolds number and geometries. Further work is on the way to extend the 
comparison to elliptic flows, Focusing the attention on those models that proved to 
possess the necessary degree of generality. 

Some of the models, like the two-layer formulations, showed a poor fit with ex- 


12 



periments. This is not surprising since these, and other, models had been tuned to 
fit high-f2e flows using different data sets. Most of these formulations proved to give 
satisfactory depictions of more complex flow patterns than the one investigated here. 
Nevertheless, the comparison of the models in such a simple, and ideal, flow case al- 
lowed highlighting some differences that would have probably been lost in a complex 
flow configuration. 
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